[N]pT Monte Carlo Simulations of the Cluster-Crystal-Forming Penetrable Sphere 

Model 



(N 

o 

< 

(N 



Kai Zhang^ and Patrick Charbonneau-'^'^ 

^Department of Chemistry, Duke University, Durham,, North Carolina, 27708, USA 
'^Department of Physics, Duke University, Durham, North Carolina, 27708, USA 

(Dated: February 27, 2013) 

Certain models with purely repulsive pair interactions can form cluster crystals with multiply- 
occupied lattice sites. Simulating these models' equilibrium properties is, however, quite challenging. 
Here, we develop an expanded isothermal-isobaric [N]pT ensemble that surmounts this problem by 
allowing both particle number and lattice spacing to fluctuate. It is particularly efficient at high 
T, where particle insertion is facile. Using this expanded ensemble and thermodynamic integration, 
we solve the phase diagram of a prototypical cluster-crystal former, the penetrable sphere model 
(PSM), and compare the results with earlier theoretical predictions. At high temperatures and 
densities, the equilibrium occupancy n%'^ of face-centered cubic (FCC) crystal increases linearly. At 
low temperatures, although n'^ plateaus at integer values, the crystal behavior changes continuously 
with density. The previously ambiguous crossover around T ~ 0.1 is resolved. 

PACS numbers: 64.70. K-,64.70.D-,82.30.Nr,62.20.-x 
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I. INTRODUCTION 



Motivated by the ubiquity of repulsive electromagnetic 
interactions, most models for matter include a core com- 
ponent that diverges at zero separation, which results 
in a first-order fluid-crystal transition. Except for a few 
exotic examples, such as certain bosonic systems [l|, y|j 
simple fluids of atoms and small molecules thus simi- 
larly crystallize. Yet when nano- and micro-scale par- 
ticles are concerned, more complex models that include 
bounded repulsive interactions are also conceivable. Soft 
and mostly empty particles, such as dendrimers, for in- 
stance, can sit on top of each other with only a finite 
energy penalty Q. A soft-core interaction captures the 
effective coarse-grained nature of the potential of mean 
force between them Q, and at relatively low densities 
these interactions are reasonably pairwise additive [5|] . 

Liquid state theory suggests that the phase behavior of 
such soft-core repulsive interactions falls into one of two 
categories, depending on the presence or not of negative 
components in the Fourier transform of the pair poten- 
tial ~ or, roughly speaking, on its steepness [y, 0|- On 
the one hand, fiat potentials with purely positive Fourier 
transforms undergo reentrant melting. Under compres- 
sion at low temperatures, the system first crystallizes and 
then remelts. The Gaussian core model (GCM), which 
was first proposed as a model for plastic crystals [g| and 
has since been thoroughly studied |9l-[l2|. exhibits such 
behavior. On the other hand, steep potentials with neg- 
ative Fourier components cluster up at high densities. 
Under compression, these systems form crystals in which 
each lattice site may be occupied by multiple particles. 
The transition from fluid to cluster crystals is driven 
not only by gains in free volume, but also by reducing 
the overlap energy [7|. The penetrable sphere model 
(PSM), which was initially proposed to describe the in- 
teraction between polymeric micelles [1^ , belongs to this 



second category. Its properties have been explored by 
density functional theory (DFT), cell theory, kinetic the- 
ory, and basic simulations [l4| - [l8| . To complement these 
two systems, Mladek et al. have introduced the gener- 
alized exponential model of index n (GEM-n), with a 
pair potential that scales with distance r like exp (— r"), 
to interpolate between the GCM (n = 2) and the PSM 
(n — cx)) 19] . Since then, a number of additional dy- 
namic and thermodynamic studies have looked at the 
exotic properties of this class of models [3, [20l - t25| . 



Studying the thermodynamics of cluster crystals by 
simulations is, however, particularly challenging. When 
initializing the system, there is a trade-off between con- 
figurations with more lattice sites but fewer overlaps, and 
configurations with fewer lattice sites but more overlaps. 
The problem is that once initialized, (reasonably) finite 
systems cannot find the average lattice site occupancy 
n°'^ because no simple ensemble allows for the lattice 
spacing and occupancy to simultaneously relax |7l. l2ll . l26j . 
Incommensurability and free energy barriers get in the 
way. A free energy scheme based on an extended ther- 
modynamic description is necessary to locate the equilib- 
rium phase, but the approaches used thus far have come 
at a fairly high computational cost [3, [2l|, [2J|. In this 
paper, we propose a more direct and efficient extended 
ensemble that effectively allows both the particle number 
N and the box volume V to fluctuate. We use it to solve 
the equilibrium properties of the canonical PSM. 



The plan of this paper is to flrst revisit the PSM and 
the basic thermodynamics of cluster crystals (Sect. HI)) 
and to introduce the simulation methodology (Sect. Hill . 
After discussing the results (Sect. HV|) . a brief conclusion 
follows. 



II. MODEL AND THERMODYNAMICS 

The penetrable sphere model (PSM) was first mtro- 
duced by Marquest and Witten as a prototype for soft 
matter interactions [Tj, [23|- It is described by a pair 
potential 



(r) = 



e for < r < (T 
for r > a, 



(1) 



where the particle diameter a sets the unit of length and 
the barrier height e > sets the unit of temperature T 
(/3 = l/T) with Boltzmann's constant ks set to unity for 
notational convenience. The model can also be seen as 
the n — )■ oo limit of the GEM-n, with pair potential 



= ..-('■/-)" 



(r) = ee 



(2) 



The phase behavior of the PSM has not been accu- 
rately determined, but some of its general features are 
reasonably well understood. Based on liquid state the- 
ory arguments, one expects that systems with bounded 
yet harshly repulsive pair interactions, such as the GEM- 
n with n > 2 [^ [l3|, should form face-centered cubic 
(FCC) cluster crystals at high T and number density 
p = N/V. Its Nc lattice sites then have an average occu- 
pancy Tic — N/Nc > 1. Density functional theory QDFT) 
predicts that ric increases linearly with p at high T [6lll7l|. 
while at low T the fluid-solid transition approaches the 
hard sphere limit [l^. A cell theory treatment further 
suggests that increasing p in the low T regime results in 
a continuous evolution of the lattice occupancy [Ij, [l^ , 
in contrast to the series of first-order isostructural transi- 
tions and critical points observed in the GEM-4 [23, l2J] . 



A. Expanded Thermodynamics 

In order to study the behavior of cluster-crystal form- 
ers, we first review and expand a thermodynamic frame- 
work for treating multiple-occupancy lattices. The core 
problem with simulating these crystals is that even for 
very large systems the average lattice site occupancy 
cannot generally relax to its equilibrium value once the 
crystal forms a lattice commensurate with the simulation 
box. It is only possible to introduce or withdraw planes 
of defects without causing elastic stress, which severely 
restricts the extent to which ric can be tuned. Inspired by 
the generalized thermodynamics of systems with vacan- 
cies and interstitials [20|, phase rules that take into ac- 
count the possibility that crystals have multiple particles 
occupying each lattice site have been formulated [3, [21| . 
In order to specify an equilibrium phase, Mladek et al. in- 
troduced an additional pair of thermodynamic variables: 
A^c and an associated chemical potential-like quantity pc, 
roughly corresponding to the free energy penalty of in- 
serting a lattice site. 

In this expanded formulation, the constrained 
Helmholtz free energy Fc{N,V,T,Nc) has an exact dif- 



ferential form 

dFc = -SdT - pdV + pdN + p^dN^, (3) 

where the chemical potential p , the pressure p, and the 
entropy S are all implicit functions of ric, and at equi- 
librium F{N,V,T) = Fc{N,V,T,N^'i). The constrained 
Gibbs free energy Gc{N,p,T,Nc) ^ pN + p^N^ has a 
similar differential form 



dGc = -SdT + Vdp + pdN + p^dN^ 



(4) 



and G{N,p, T) = Gc(iV,p, T, N^'i). Because p and p^ are 
functions of A^c, they do not correspond to the system's 
equilibrium quantities unless iVc minimizes the overall 
free energy of the system, i.e.. 



fdFc 



\dNc 



fdG, 



c/ Ny,T;N^=N':'^ 



\dNc 



c / N,p,T;N^=NS" 



= 0; or ( TTTT ) = 0- 

(5) 

Equilibrium thus occurs when the driving force for chang- 
ing the number lattice sites vanishes, i.e., pc ~ 0. It 
is, however, much more computationally convenient to 
fix Nc and vary N. The free energy minimization is 
then applied on the per-particle quantities, /c(p, T, Uc) = 
Fc{N,V,T,Nc)/N or gc{p,T,nc) = Gc{N,p,T,Nc)/N, 
such that 



dfc 
duc 



0; or . 



dgc 



0, (6) 



as obtained using a generalized Gibbs-Duhem relation. 

B. Response Function 

This thermodynamic construction provides important 
physical insights into the equilibrium response functions 
of cluster crystals. The heat capacity Cy, bulk modulus 



(or compressibility) B( 



) and coefficient of thermal 



expansion a, for instance, can be described as having two 
different physical contributions. The first comes from 
"normal" crystal relaxation, i.e., an affine transforma- 
tion under compression, and the second from clustering. 
For reasons that will become clearer below, the affine 
component is also dubbed the virial contribution. 

The response function is generally defined as the 
derivative of a thermal quantity Y (energy, enthalpy, vol- 
ume, etc.) with respect to a changing parameter x (tem- 
perature, pressure, chemical potential, etc.). In cluster 
crystals, Y itself is a function of N^'^, which in turn de- 
pends on X, so 



dY{x, N^'i) 
dx 



dY 
dx 



N^=N°'^ 



dY 



x:N^=N, 



(dNl\ 
'-ci\ dx 1 



where the two terms correspond to the virial and cluster- 
ing mechanisms, respectively J28l|. The bulk modulus of 



cluster crystals, for instance, 
^^ ^/ dp{N,V,T,N^^) \ 



V 



V dV 

dp 



ON, 



/ N.T 



= -V 



dp 
dV 



N.T,N^=N''^ 



c/ N,V,T-N^=N^'^ V C'V' 



(7) 



N,T 



has an additional softening contribution due to cluster- 
ing. Under compression, the lattice relaxes not only by 
affinely reducing the lattice constant, but also by elim- 
inating lattice sites and clustering [3, ^M, [2J|- Analo- 
gously, the coefficient of thermal expansion 



f fdV{N,p,T,N^'i) 



V 



1 r dv 



dT 



N.p 



V [df 



N,p.N^ = N°'^ 



V \dNc 



c / N,p,T;N^=N^ 



(dN^\ 



(8) 



N,p 



reveals that under isobaric heating the dilation of the lat- 
tice is accompanied by a break up of the clusters them- 
selves, which consumes some of the additional expansion. 
Instead of calculating the thermal derivative, it is well 
known that response functions can also be expressed as 
fluctuations of thermal quantities, e.g., Cy = (3{{E'^) — 
{E)^)/T 29, 30]. Ahhough formally also true for cluster 
crystals, special attention needs to be given to Nc being 
held fixed in simulations. In that case, the calculated 
fluctuations only capture the affine contribution from the 
lattice. In other words, the simulation at fixed N^ does 
not measure the fluctuation of equilibrium lattices, but of 
lattices constrained to the specific equilibrium occupancy 
at one phase point. The clustering contribution must 
thus be separately calculated, in order to calculate the 
correct equilibrium behavior. 



III. METHODS 

Traditional simulation methods either fix the number 
of particles, as in constant NVT or NpT simulations, 
or the lattice spacing, as in constant fxVT simulations, 
which in all cases prevent cluster crystals from relaxing to 
their equilibrium ric- An ensemble that would allow both 
N and V to fluctuate may naively seem to resolve the 
problem, but such /ipT ensemble would have unbounded 
fluctuations. Earlier simulation approaches thus em- 
ployed elaborate free energy schemes in order to minimize 
/c and (7c ,24] or equivalently to locate /Xc = [3, [2l| . In 
these schemes, locating the equilibrium structure requires 
tens of free energy calculations, each necessitating one to 
two dozens of independently sampled state points. This 
high computational cost limits the method's applicabil- 
ity to all but the simplest models. In order to reduce 
the computational burden and broaden the range of ac- 
cessible models one should limit resorting to free energy 



integrations. We present below such an approach based 
on reformulating the expanded isothermal-isobaric [N]pT 
ensemble ISll. 



A. [N]pT Ensemble 

To elucidate the [N]pT ensemble, it is useful to flrst 
review the generalized (great grand canonical) ensem- 
ble [S^, |3J] . Its trivial partition function 



T 



N 



V E 



pE 



(9) 



N 



yy^^^{N,p,T)^ 



N 



V 



^{^i,v,T) = l, 



derived here using a discrete state (histogram) formal- 
ism even for continuous variables, implicitly defines the 
microcanonical fi, canonical Q, isothermal-isobaric A, 
and grand canonical S partition functions. As mentioned 
above, in this ^pT ensemble no extensive quantity is spec- 
ified, which results in N ^V and E having unbounded fiuc- 
tuations. 

An [N]pT ensemble can solve this last problem by con- 
straining N within bounds [A'min, A^max] that (if properly 
chosen) enclose rf^ at the studied p, T, and Nc- The new 
partition function 



'Y^* 



A'max 

>: 


g/3G^ 


V 


-fipV \ ^ 
E 



n{N,V,E)e- 



-I3E 



(10) 



has weights Gn = g^N that are both explicit (extensive) 
and implicit functions of TV to control the probability of 
visiting states of varying N . These weights are not known 
a priori, but can be self-consistently determined, as is 
done for the multicanonical ensemble j34j , the expanded 
ensemble parallel tempering [3^] , and Wang-Landau sam- 
pling i3a, [33 • Once the weights are determined, the joint 
probability of observing a state with a specific N, V, and 
£;is 

nN, V, E) = l^n{N, F, ij)e/5GN-/3py-0i=;^ ^^^^ 

and the marginal probabilities of observing a state with 
specific N and V or with a specific N only are, respec- 
tively. 



^ V{N, V, E) ~ e'^^"-''P^Q(iV, V, T) 

E 

^N(A^)-^^^(A^,T^,i?)^e^^-A(iV,p,r). (12) 



E 



E V 



1. Simulation Method 

During a simulation, in addition to standard particle 
and logarithmic volume displacements [38] , particle inser- 



tions (+) and removals (— ) are used. In order to preserve 
detailed balance, the acceptance ratios are 

acc(7V, V,E^E + AE) = min {l, e'^^^} (13) 

acc(iV, V^V + AV,E^E + AE) = 
min { 1, e-^(A^+pA^)+(^+i) '" ^^^^ } (14) 

acc(iV ^ N±1,V,E-^ E + AE^) = 
min{l,,7±e/'^«*-^^^*}, (15) 

where AG± = Gn±i - Gn, r/+ = V/{N + 1), r]- ^ 
N/V, and AE^ is the energy cost of inserting/removing 
a particle. 

In order to self-consistently determine Gat, for each 
iteration i we obtain a histogram W* (N, V, E) that keeps 
track of the frequency at which each [N, V, E) state is 
observed. After i iterations, the normalized histogram 
gives the probability of observing a state 



V'iN,V,E) 



W'{N,V,E) 



(16) 



EN,V,E^'iN,V,E) 

that approximates the density of state 

n\N, V, E) = V'{N, V, ^)e-'3G^+/3pV+,3£; ^^^^ 

and the other partition functions 

Q'{N, V,T)^Y^ n'{N, V, E)e~'^'^ and (18) 

E 

A%N,p,T) = J2Q\N,V,T)e-^P''. (19) 

V 

Note, however, that even once the scheme has converged 
n, Q and A all lack a multiplicative constant T*, so the 
conjugate field that is used for the (z + l)th iteration 



/3GJ+1 



lnA*(7V,p,T) 



(20) 



also lacks an additive constant G*. The iterative proce- 
dure is repeated until the marginal distribution Vn (N) is 

flat within [iVminjiVmax], i.e., ef'^^ A{N,p,T) = constant 
(Eq. [T^ , which coincides with the convergence criterion 
Gn^' = G^ + G' for aU N e [iV^i„, N,^,^]. 

A final simulation using the converged weights con- 
structs the equilibrium W{N, V, E). This last [N]pT sim- 
ulation is conceptually equivalent to running a series of 
NpT simulations with different iV's, where the transition 
probability of hopping from one to the other is governed 
by Gn- The weights thus capture the constrained Gibbs 
free energy Gc of the system up to a constant G. Be- 
cause the equilibrium occupancy is determined by mini- 
mizing Qc with respect to N (Eq. |6]) , the unknown con- 
stant G, which leads to a C jN term, must be obtained by 
thermodynamic integration at a given A'o G [A^min, A^max] 
(Sect.|IIlB|. 
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FIG. 1: [N]pT results for the FCC cluster crystal with Ac = 
500 at r = 1.5 and p = 33.4. The weights are initialized 
with the chemical potential p at Ao = 3910 (dashed line) 
and converge to ^n = 0.00078065A -I- 9.6940 -I- 11877/7V (0). 
The constrained Gibbs free energy per particle (solid line), 
after adding the constant shift C — 397, is minimal at N"'^ = 
3965(10). Inset: The converged distribution Vn is flat within 
the sampled bounds, fluctuating around 1/151 (black line). 



2. Simulation Details 

For high temperature simulations, 4 x 10"* Monte Carlo 
(MC) cycles are performed for A^ = 1000-5000 particles 
at each iteration. Each cycle contains on average 30%A^ 
particle number changes, 70%N particle displacements 
and one volume fluctuation. For the crystal phase, the 
particles are spread over A^c = 256-1372 lattice sites de- 
pending on the average occupancy. A particle range of 
AA^ = A^max - A^min = 150~200 is found to provide a 
reasonable trade off between the need to correctly esti- 
mate the equilibrium occupancy and the computational 
cost. The histogram of states separates volume fluctu- 
ations in 50 bins and energy fluctuations in 500 bins 
that span the observed range of fluctuations. For most 
system sizes studied, the range has Knax — Knin ^ 50 



and En 



E„ 



2000. The summation over these 



states is done on logarithmic scale to avoid numerical 
overflow |39| . 

A good initial guess for the weights is the chemical po- 
tential of a system with A^q S [A^min, A^max], 5n = m(^o), 
obtained by Widom insertion in a short NpT simula- 
tion [38[. From this point on, fewer than ten iterations 
are needed to converge the standard deviation of the 
equilibrium probability distribution within cr-p„ ~ 0.001. 
Even though the resulting weights are noisy for a given 
A^, considering the full A^ range, by, for instance, fitting 
them with a parabola smooths out these uncorrelated 
fluctuations. After including the offset G, the constrained 
Gibbs free energy thus has an approximate form 

Gc -- C2N^ + ciN + Co, 

where co,ci, and C2 are numerically determined. Corre- 



spondingly, 

ffc = C2N + ci + co/N, 
^ — 2C2N + ci, and 
C2 ,.2 , Co 

If appropriate bounds are chosen, g^ is then minimized at 
N""^ = Vco/c2 e [A^min,A^max]- Figure [H iUustratcs this 
operation. Starting with fj,{No = 3910) = 15.784672, 
the simulation converges to ^n- After including the con- 
stant offset C = 397, gc is minimal at N'^'i = 3965(10) 
{ric = 7.93(1)), where the error is obtained from running 
several optimizations of 5n- The equilibrium particle 
number obtained by this method is in good agreement 
with that determined by pure thermodynamic integra- 
tion, N^'^ = 3960(10), following the approach described 
in Ref. 24] and Section IIIIBI This particular example 
also illustrates that the minimum of gc does not need 
correspond to that of 5ni and can even be out of the cho- 
sen range [A'^jn, iVmax], although for the minimization to 
be reliable within the accuracy reported here the range 
should be chosen such that the minimum does fall within 
those bounds. 

Determining the equilibrium behavior from the free en- 
ergy requires a fairly high degree of numerical accuracy. 
It is thus worth estimating the error made in this process. 
The typical change in weight values over the full N inter- 
val A(7c, is smaller than the shift from ^n to gc. For in- 
stance, over an interval of length AA^ ^ 10^, Age ^ 10^"^. 
The Gibbs free energy per particle, which is of order 10^, 
is determined by thermodynamic integration with an er- 
ror 6 go ^ 10~^. The offset constant C will thus have an 
error SC = SgoNo ^ 10"^ x 10^ = 10^ - Scq. The error 
Sgc in g^ is then of order 10^^, which is an order of mag- 
nitude larger than the variation A^c- Yet the position of 
the minimum of the curve, iV"'^ — -\/co/c2 is only shifted 
by SN'^'i - N°%^1 + Sco/co - 1) ~ 10^ x lO'^ = lO". 
Therefore, although the error in the vertical shift of gc 
may be an order of magnitude larger than the shallow- 
ness of the curve itself, the position of the minimum is 
preserved with an accuracy comparable to that obtained 
from running multiple instances of the optimization. 



3. Low-Temperature Implementation 

Special attention needs to be drawn to use the [N]pT 
method for the PSM at relatively low temperatures, 
T < 0.1. First, the traditional volume fluctuation al- 
gorithm, in which the positions of all the particles are 
affinely rescaled, is inefficient to simulate the PSM. Low 
temperature compressions of the PSM, like in a system of 
dense hard spheres, result in additional particle overlaps, 
whose cost is prohibitively high. To overcome this prob- 
lem, we use an alternative volume fluctuation algorithm 
recently developed by Schultz and Kofke [40|, which im- 
proves the sampling efficiency by nearly three orders of 



magnitude at T = 0.05. In the Schultz-Kofke algorithm, 
when volume changes from V to V^ -I- AV^ on logarithmic 
scale, the position of each lattice site Ri is first scaled 
affinely as Rj,„ow = R.^,old[{V + AV)/V]^^^ , then the 
distance of each particle from its lattice site Vj is changed 



V 



V + AV 



l/{3N-3) 
,/3(pAV+AEiat) \ ^ /2I) 



The modified acceptance ratio is then 

acc(y ^ V+AV,E^ E+AE) = min |l,e"'^('^-^"'^^""H , 

(22) 
in which the effect of volume change is redistributed into 
the Tj scaling step [42, and where ASiat is the change 
in lattice energy in the Ri scaling step. For step-wise 
discontinuous potentials like the PSM, Aii^iat is precisely 
zero if the particles sitting perfectly on lattice sites do 
not overlap before nor after a volume displacement. Yet 
if Ai^iat were nonzero, it would scale like A^, because all 
the particles on neighboring lattices would go from not 
overlapping to overlapping. Volume compressions with 
Ai?iat > are thus rejected with overwhelmingly high 
probability at low temperatures. 

Particle insertion also becomes difficult at low temper- 
atures. For low T crystals with an occupancy interme- 
diate between two integers, new particles can only be 
successfully inserted at the lower-occupancy positions. 
Modified particle insertion schemes, such as the staged- 
insertion algorithm designed for fiuid phase [42, |4j| , do 
not help, because the problem is intrinsic to the system. 
Inserting particles requires one to randomly identify a 
vacancy in the crystal, which constrains the [N]pT ap- 
proach at low T to crystals with non-integer occupancy. 
Once these regions are identified, PSM crystals with inte- 
ger occupancy can then be treated using standard NpT 
simulations, for instance. For the [N]pT scheme at low T, 
the chemical potential used as initial guess for Gn is ob- 
tained from staged Widom insertion [42, 44], and 10^ MC 
cycles are needed at each iteration to obtain good statis- 
tics. To speed up the simulations, more particle inser- 
tions and deletions (50%A^) and ten to a hundred times 
fewer volume fluctuations are used in each cycle. The 
particle number window is also shortened to AN = 50. 
Yet the inherently low efficiency of particle insertion in 
low T cluster crystals leaves the [N]pT approach less effi- 
cient than the free energy integration scheme of Ref. [2J. 
Only a few state points were thus equilibrated this way. 
The rest were obtained with the traditional approach 
that involves multiple thermodynamic integrations [24i] . 



B. Thermodynamic Integration 

Free energy calculation from Kirkwood thermody- 
namic integration couples the system of interest to a 




enough fluctuations and a smooth integration path. The 
constrained free energy of the sohd phase with a certain 
occupancy is then 



/c = /o + 



($ 



*0>A 



N 



dX, 



(28) 



as illustrated in Fig. [2] 46]. It is essential for the inte- 
gration path to be continuous and tp show no hysteresis, 
as the preservation of crystal symmetry in Fig. [2] also il- 
lustrates. For fixed p,T,Nc, free energy integrations are 
performed for various iV's until the equilibrium nl'^ is 
identified from Eq. [S] [2J] . 

The free energy of the fluid is obtained from thermody- 
namic integration from infinite temperature limit where 
the system is effectively an ideal gas [38| 



FIG. 2: Top: Therniodynamic integration of the FCC solid 
phase with Wc = 2580/1372 = 1.88 at T = 0.1 and 
p = 1.8. Bottom: The radial distribution function g{r) = 
~(X]i=2 ^i^~^i)) ^t various A shows that the crystal order is 
preserved along the integration path. 



reference system whose free energy can be obtained ex- 
actly [3a, |43- For the cluster crystal, the potential 
energy of the coupled system with particles at r^ = 
{ri,r2,- • ^rjv} is defined as in Ref. vA \4B 



^N\ 



<^x^{l-\)Mr'l + mr'l 



where 



N~l N 



i—1 j>i 



(23) 



(24) 



and the coupling parameter A varies from to 1. We 
choose the reference system to be free particles under 
spherically attractive potential wells centered at the lat- 
tice sites 



with 



Mr) = 



Mr'l 



_ J eo < 




N 

E 



M^i), 



re UiI\«o(R-j 
otherwise 



(25) 



(26) 



where wo(Ri) is the size of one attractive spherical well at 
lattice vector R^ . The free energy of the reference system 
at temperature T is therefore [46[ 



/3/o = In 



N 



(V - Vo) + Voe^f^^o 



1, 



(27) 



where Vq — NcVq is the total volume of the pote ntial 
wells and the thermal wavelength A = yT^/i^/^Trm is a 
constant at fixed T and is here set to unity. The param- 
eters Vq and eo are tuned for each simulation to obtain 



Pf^pr+ / {'^)/Ndp' 



(29) 



which is equivalent to setting eo = in the crystal inte- 
gration scheme. Standard integration over density is also 
used 



/(P2) = /(Pl) + / ^dp. 
J Pi P 



P2 



P 



(30) 



to verify the results. 



C. Histogram Reweighting 

Because [N]pT ensemble simulations rely on building 
the full density of state, histogram reweighting provides 
the thermodynamic properties of the system at neighbor- 
ing T andp [3|. 



Reweighting Over Pressure 



Histogram reweighting over pressure, for instance, pro- 
vides -Bvir and I -5^ ) from a single [N]pT simulation. 

For each N within the bounds, the volume V is sampled 
with the conditional probability, knowing N, at pressure 
P 



rv\N{v\N) = 



Vnv{N,V) _ e'f'P^Q{N,V,T) 



Vn{N) 



A{N,p,T) 



(31) 



The thermal average of M configurations obtained from 
Monte Carlo sampling at pressure p is thus 



(V) 






M 



(32) 



Now, suppose the simulation is run at a nearby different 
pressure po- The probability of observing a state becomes 



Vv\Nfi{y\N) 



e-I^P«^Q{N,V,T) 
^{N,Po,T) ■ 



(33) 



In order to evaluate the same average quantity (V) at p, 
while sampling configurations at Pq subject to the distri- 
bution VviN.Oi the histogram is reweighted 



{V) = 



EvVVylj^^^e-^P^'Q _ ^^ye-/5(J'-P«)^ 



E 



vT'vko''-^'"'Q 



E. 



-f3{p-po)V 



(34) 



For each Uc, i.e., for each N at fixed TVc, this re weight- 
ing scheme provides (v) = {V)/N over a short pressure 
interval enclosing pq. By integrating (v) for each ric, the 
constrained Gibbs free energy per particle over a pressure 
interval follows 



5c(p,"-c) - gc{PQ,nc 



{v)dp'. 



(35) 



For each neighboring pressure p, gc{p,nc) is minimal 
at the equilibrium occupancy n'^'^ and the equilibrium 
Gibbs free energy per particle g{p). The corresponding 
equilibrium Helmholtz free energy per particle fdp) = 
9c{p) — p{v) is similarly available for neighboring densi- 
ties. Numerically, the volume for each N obtained from 
[N]pT simulations is reweighted at neighboring pressures 
and then fitted to a third-order polynomial in pressure, 
before using Eq.[35l Figure [3] shows an example of the re- 
sulting 5c(p, n-c). As a rule of thumb, histogram reweight- 
ing is sufficiently accurate as long as N°^ remains within 

the bounds [A^min, A^max]. 



2. Reweighting Over Temperature 

Histogram reweighting over temperature similarly pro- 
vides ttvir and ( -^f- ] from a single [N]pT simulation. 

One can show that the thermal average of the volume 
at neighboring temperatures can be calculated from the 
configurations sampled at Tq 



(V) 



J2v Ve-^I^-Mpv j2j^ n{N, V, E)e-^l^-l^°^^ 



E^ 



-(13- 



■MpV ^^ n{N, V, E)e-(I^-ME 



(36) 

To calculate the Gibbs free energy at fixed p and neigh- 
boring T, one reweights the isobaric-isothermal partition 
function. Although the real partition function differs 
from the sampling result by a multiplicative constant C* , 
this constant is determined by the free energy at Tq by 

PoG,{N,p,To) = -i\nAiN,p,To)+lnC*). (37) 

In a Monte Carlo scheme, the fraction phase space sam- 
pled A{N,p,T()) is proportional to the number of ac- 
cepted configurations Af, because these configurations 
are already Boltzmann weighted. Note that in this for- 
mulation, C* absorbs the normalization by the total 
number of the samples considered. For instance, if the 
algorithm is ran twice as long, C* is halved. The sam- 
pling partition function at neighboring temperatures is 
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FIG. 3: Histogram reweighting g^ at (top) T = 1.5, around 
po = 32.4, and at (bottom) p = 206 around To = 4. A 
constant /S.g has been subtracted from the free energy in order 
to superimpose the curves at the beginning of the interval. 



thus 

V E 

(38) 
and the corresponding constrained Gibbs free energy for 
each N is 

/3Gc(A^,P,r) = -(lnA(iV,p,T)-|-lnC*). (39) 

Sample reweighting results are shown in Fig. |31 

D. Phase coexistence 

Identifying the phase coexistence region requires the 
free energies of the cluster crystal and fluid phases nearby. 
For the solid phase, the free energy is obtained at various 
temperatures using [N]pT simulations, thermodynamic 
integration, and histogram reweighting, as described in 
the previous sections. The free energy of the fluid phase, 
which does not suffer from the same occupancy ambigu- 
ity, is straightforwardly obtained from standard free en- 



ergy integration techniques 38^. For 0.07 < T < 0.1 near 
coexistence, however, standard Monte Carlo sampling of 
the fluid phase results in long correlations between the 
configurations. The few particle overlaps present in those 
conditions do not easily relax because the dense fluid only 
rarely opens large enough gaps. In this regime, longer 
simulations with 10^ MC cycles are necessary to estimate 
the equilibrium quantities. 

The fluid-solid phase coexistence at a specified T is 
obtained by common tangent construction of the trans- 
formed free energy f — fp — Kp. By an appropriate 
choice of constant K, the curvature of the Helmholtz free 
energy at coexistence can be enhanced [2l|. Because 

a common tangent construction also provides the chem- 
ical potential and pressure of the two phases concerned. 
In order to better trace the coexistence line on the phase 
diagrams, its slope is also calculated at a few points. For 
the p-T projection, the standard Clausius-Clapeyron re- 
lation can be used 



dp 
df 



As 



(40) 



where Aw and As are, respectively, the per particle 
volume and entropy difference between the two phases 
across the transition. The latter quantity is obtained 
from the free energy results. For the T-p projection, each 
phase has 



dp 
df 



dp 
-pa + Pi^[-Qf 



(41) 



For the fluid branch, the response functions are directly 
obtained from fluctuation identities [30]. The compress- 
ibility is accessible from fluctuations in N for pVT sim- 
ulations, or in V for NpT simulations 



1 
V 



dV\ 
dP L 



fiV 



(iV2) - (7V)2 



PiV) 



(7V)2 



{vy 



{V)' 



(42) 



while the coefficient of thermal expansion is estimated 
from the covariance between V and enthalpy H = E+pV 
for NpT simulations 



1 
a = -— 



dV\ 
V \'df . 



nVH)-{V){H) 



T 



{V) 



(43) 



As discussed in Sec. IIIBl for cluster crystals the fluctua- 
tion (virial) contribution is insufficient. The full deriva- 
tives for the solid branch is instead obtained by nu- 
merically differentiating the equilibrium results or by 
reweighting the [N]pT simulation results performed at 
the equilibrium pressure and temperature as discussed 
in Sec. Unci 



IV. RESULTS AND DISCUSSION 

Using the methods described in the previous section 
the structure and properties of the crystal phase and of 
the fluid-crystal coexistence of the PSM model are exam- 
ined. 



A. Equilibrium Cluster Crystal Properties 

We flrst identify the equilibrium cluster crystal occu- 
pancy n°^ at different phase points. Qualitatively, n'^^ 
increases with density and pressure, and at flxed pres- 
sure the number of particles per cluster decreases with 
increasing temperature as new lattice sites are formed. 
Quantitatively, at flxed high T, n"'^ grows roughly lin- 
early with density, while at fixed density, it grows linearly 
with temperature (Fig|4l[a)). If viewed as a function of 
p and T, n^''(/9, T) has a bilinear form n'^ = aT -\-hp -\- c 
with positive a and 6, and a small c. The DFT predic- 
tion, n°^ = 2T + p [l3| , is qualitatively comparable with 
the coefficients obtained from fitting the simulation re- 
sults (FigSJa)), but in order to capture the clustering 
behavior for temperatures as low as T ^ 0.1, quadratic 
corrections are necessary. 

As T is further lowered, n°'^ plateaus at integer oc- 
cupancy values, which we denote FCC(nc). At T = 0, 
however FCC(nc) with integer ric > 2 is again only sta- 



ble at a single density p ~ 



UcP 



cpi 



where 



-'cp 



= V2 



IS 



the crystal close packing density of hard spheres. For 
ncPcp < P < [nc + l)Pcp the ground state is a random 
mixture of sites with occupancy ric and rzc -|- 1, because 
all cluster arrangements are degenerate (FiglUJb)). 

What about intermediate temperatures? In the GEM- 
4, transitions between cluster phases with different inte- 
ger occupancy FCC(7T,c)-FCC(nc -t- 1) was found to be- 
come first order below an isotructural critical point j24| . 
Although a simple phonon analysis suggests that a sim- 
ilar type of transition may be possible in any cluster- 
forming system at low enough temperature [23|, no such 
transition is observed here. This point deserves further 
consideration. For the small integer occupancy studied, a 
multiply-occupied lattice site is effectively a large coarse- 
grained particle whose size is related to the lattice con- 
stant a. To obtain a first order isotructural transition, 
i.e., demixing of FCC(nc) and FCC(nc -I- 1) phases, a 
large enough size heterogeneity between sites of differ- 
ent occupancy is necessary for the free volume gained to 
overcome the entropy of mixing, as in binary hard sphere 
mixtures [47| . Otherwise, the sites with occupancy ric 
and Tic + l randomly mix. For the PSM, the evolution of 
the lattice constant in Fig.|4l^b) indicates that the largest 
lattice difference between FCC(l) and FCC(2) is smaller 
than 3% at T = 0.05. Because this difference goes to zero 
at r = 0, its relatively rapid shrinking with T may be 
sufficient to prevent isostructural phase separations, al- 
though it cannot be excluded altogether. This qualitative 
distinction between the GEM-4 and the PSM reflects the 
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FIG. 4: (a): Equilibrium occupancy at T > 3 has a bilinear form n°'^{p,T) = 1.64(3)r + 1.132(12)p - 0.087(108). When the 
lower temperature results are included, higher order terms are necessary to obtain a reasonable fit n^'^{p,T) = 2.65(ll)r — 
0.12(1)t2 + 0.77(5)/9 + 0.014(2)p2 +0.13(9) (right inset) and nl'i{p,T) = -0.000057(7)^^ + 0.083(4)p + 8.9(4) - 0.97(ll)r (left 
inset), (b): At low temperatures, the equilibrium occupancy plateaus at integer values (bottom), but the difference in lattice 
constant a = (\/2»t-c/p)^ between neighboring integer occupancy phases is too small to induce demixing (top). At T = 0, the 
lattice constant for p > pcp is fixed and the occupancy increases linearly, (c): Relative particle number fluctuation for various 
temperatures. The theoretical prediction for the fluctuations between [njr'^j and \n'^'^] (solid line - see text) agree with the low 
temperature results. 



sensitive dependence of isostructural transitions on the 
details and convexity of the interaction potential [4^. 
The relatively long tail of the GEM-4 interaction indeed 
results in an additional energetic contribution to the lat- 
tice constant, and overcomes the entropy of mixing at 
low enough T. For GEM-n with 4 < n < oo, one ex- 
pects the threshold lattice constant difference between 
lattices with occupancy ric and ric -I- 1 to also depend on 
T via the entropy of mixing. Qualitatively, we expect the 
isostructural critical point to steadily decrease with the 
increase of n and ric. A better knowledge of the drive 
for binary hard sphere mixtures of different size ratios 
to demix would be needed to estimate the equivalent be- 
havior in the GEM-n. 

With the equilibrium occupancy at hand, we study 
how the number of particles on a given site fluctuates in 
various T regimes. This quantity is particularly impor- 
tant in determining which of the basic theoretical approx- 
imation is most analytically reasonable. Overall the rel- 
ative fluctuation decreases with increasing temperature 
and density, but this progression is not monotonic (Fig. |4] 
Right). At low temperatures (T < 0.1), a system with an 
equilibrium occupancy intermediate between two integers 
[n""^] < 71°'^ < [n""^] can be seen as a random mixture of 
sites with occupancy [n'^^\ and occupancy \n^'^'] . Exci- 
tations creating \nl'^~\ -|- 1 or [n^'^\ — 1 defects are strongly 
suppressed. In that limit, the problem can be solved an- 
alytically. The probability of a site having [n°'^J ( [n°'^] ) 
particles is n'^'^ — [nji'^j ([ric'^] — "-c'^), and the resulting 
variance is al^ = (SL^^^J + IX"* - L^'^J r^"*! " «'*)^- 
The T = 0.05 results agree very well with this approxima- 
tion (Fig. lU^c)). This behavior explains why the mean- 
field cell model works increasingly well with decreasing 
T and nearly quantitatively for T < 0.03. As temper- 
ature increases, fluctuations become more pronounced, 
and vary continuously with occupancy. In this limit the 
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FIG. 5: The bulk modulus B and its virial contribution Bvir 
(Eq.gU at r = 4 and T = 0.05 (inset). The virial part B^ir 
is significantly softened by the cluster contribution, except in 
the low T plateau regime (p — 2.4). 



crystal behavior is better approximated by a DFT-type 
treatment in which the free energy is optimized with re- 
spect to occupancy and fluctuations, assuming a Gaus- 
sian distribution [17| . 

As determined in earlier studies [2l|, [2J| , the bulk mod- 
ulus of cluster crystals, once decomposed into its two con- 
stituents (Sect. HTB)) . shows a strong softening correction 
due to the creation/annihilation process of lattice sites at 
high temperatures (Fig. [5]) . But in the integer occupancy 
regime, such as at T = 0.05 for p — 2.4, the virial part 
is the full resistance to external compression [2J|. The 
two contributions to the coefficient of thermal expansion 
can be similarly calculated. At T = 4 and p = 206, 
for instance, calculations of the different components of 
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FIG. 6: Left: Simulation T-p phase diagram obtained from the coexistence points (O) and slopes (short black solid line). The 
other solid lines are guides for the eye, and the fluid-solid coexistence is shaded in grey. The errors are smaller than the symbol 
size. The low T scenarios offered by the cell model [14l | (dashed line) and DFT [I^l (dotted line) are shown for comparison. The 
boundaries of integer occupancy crystal phases (red solid line) terminate at T = results (•) and compare nicely with the cell 
model for T < 0.03, but not above. The * region is discussed in the text. The agreement of the high T results with DFT [13] 
(A) break down for T < 3.0 (inset). Top Right: n°'' at melting shows a rapid crossover between 0.07 < T < 0.1 that connects 
the single-occupancy regime at low T and the linear growth regime at high T (inset). Bottom Right: The corresponding 
crossover in the p-T coexistence line separates the high temperature regime (inset) and the hard sphere-like regime. In this 
last region the coexistence pressure closely follows the hard sphere value Pcoox = 11.5767 [i^] (dashed line). 



TABLE I 



Calculation of a at T 



4.0 and p = 206 with 
nr = 19.77 using (al) a fit of the T = 3.0-5.0 results (Fig.©, 
(a2) multiple NpT simulations with ric = 19.77, (a3) multiple 
NpT simulations with n^'^ determined in (al), (61) histogram 
reweighting (Fig. [3]), (62) histogram reweighting (Eq. [36]), (c) 
a direct [N]pT simulation (Eq. I32|) . and (d) the fluctuation 
identity fEa. [43)) . 



0.195(10)"^ 
0.202(10)"+''2+" 



BT 



0.169(10)"^ 
0.170(10)''^ 
0.168(10)'' 



0.97(10)"^ 

1.22(10)" -0.0265(3) 



Eq. [8] agree with each other and with the decomposition 
(Tabled). 



B. Phase Diagram 

The fluid-crystal transition, which is the only phase 
transition in the PSM, is extracted from the equilibrium 
free energies of both phases. The high and low T regimes 
are treated separately, because the underlying physics is 
markedly different above and below the onset of cluster- 
ing at coexistence around T ^ 0.1. 

At high T, the coexistence regime on the T-p phase 
diagram and the coexistence pressure on the p-T phase 
diagram decreases smoothly towards T ^ 0.1. At melting 
'^c'^(Pcocx) increases roughly linearly with T, because the 
coexistence density itself increases roughly linearly with 
T, as discussed above (Fig. [5]). DFT results capture the 
transition relatively well for T > 3, but strongly deviate 



below that point [T3| . Yet attributing the deviation to a 
stronger hard-sphere character below that T [ij], is not 
supported by the numerical results, so another physical 
assumption of the theory may then be breaking down. 

For r < 0.1 at coexistence, two conflicting descriptions 
have previously been suggested (Fig. ^. First, a mixed 
DFT (fluid) and cell theory (crystal) study predicted that 
the fluid-solid transition continuously approaches the HS 
limit upon lowering T [1J| . Although the solid free energy 
was found to be indistinguishable from the HS results at 
finite temperatures, changes to the fluid phase resulted 
in a steady drift away from the HS behavior for T > 0. 
Second, a full DFT treatment predicted a crossover tem- 
perature from a HS-like coexistence to a different regime 
around T = 0.35 [1^. We flnd the second scenario to 
be qualitatively correct, although its predicted crossover 
temperature is nearly three times larger than the numer- 
ically determined one. For T ~ 0.1, the coexistence crys- 
tal density gets smaller than the HS close packing den- 
sity Pep = V^- Clustering then gets suppressed and the 
fluid-solid coexistence curve is inflected. Below T < 0.07, 
the fluid coexists with an essentially singly-occupied FCC 
phase. The transition is hard sphere-like with the corre- 
sponding coexistence pressure and densities (Fig. |6|). 

That there should be a crossover can physically be un- 
derstood from the fact that at T ^ 1 particle overlaps 
are rare and uncorrclated. They therefore marginally sta- 
bilize both the fluid and crystal phases in a similar way. 
Once the concentration of overlaps becomes sufficiently 
high, however, the difference in structure between the 
ordered and disordered phases matters. This distinction 
between the high and low T regime also has a dynamical 
signature [la, [l^. Collisions between particles can be 
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divided in two types: soft refractive collisions, in which 
a particle goes through another, and hard reflective col- 
lisions, in which particles elastically bounce back from 
each other. At temperatures T < 0.3 from molecular 
dynamics simulations [l8| or T 5. 0-25 — 0.5 from an 
Enskog-type theoretical analysis [ia|, the first collision 
type is highly suppressed because the particle momenta 
are low, and the collision frequency of the second type is 
as high as in hard spheres. 

Unlike for the GEM-4 model, where first-order phase 
transitions between integer-occupancy lattices are ob- 
served at low T, crystal lattice occupancy in the PSM 
changes continuously. The stability regime of integer- 
occupancy phases, such as FCC2 and FCC3, only oc- 
cupies a narrow slice of the phase diagram at low tem- 
peratures. Above T ~ 0.1, these phases join the con- 
tinuum of non-integer-occupancy crystals whose T-p sta- 
bility regime is but a line that extends smoothly up to 
the liquid-solid coexistence boundary. The cell model 
fails to capture this regime's behavior (Fig. 15]). It no- 
tably erroneously predict the existence of an intermedi- 
ate single-double-triple occupancy regime (* in Fig. [S]). 
The limited extent of occupancy fluctuations at this low 
T (Fig. mjc)) shows that such high-energy excitations 
are unphysical. Simulations indeed clearly delineate the 
single-double from the double-triple occupancy regimes. 



V. CONCLUSION 

In this paper, we have reformulated the extended ther- 
modynamics of cluster crystals such that it lends it- 



self to [N]pT ensemble simulations [ST] . Using this ap- 
proach and thermodynamic integration, we efhciently de- 
termined the phase diagram of the canonical penetrable 
sphere model, which forms cluster crystals at high densi- 
ties. The [N]pT ensemble approach naturally allows for 
histogram reweighting, which was here used to calculate 
the response functions, and could be particularly useful 
for studying the critical properties of cluster crystal for- 
mers, such as the GEM-4. 

The resulting formalism and method are analogous 
to a constant pressure version of the approach used for 
determining the binary hard sphere mixture phase di- 
agram [43, [5^, which it should be able to calculate 
with some adjustments. The approach also allows for 
the study of models that form crystals in which the 
vacancy concentration is relatively large, such as hard 
cubes [5l|, [S^] ■ The key limitation for these implementa- 
tions is the need for an efficient particle insertion algo- 
rithm and a good initial guess of the conjugate field. 



Acknow^ledgments 

We gratefully thank N. Wilding for numerous dis- 
cussions and particularly for the key suggestion that 
the [N]pT approach may be useful to the study of 
cluster crystals. We also thank D. Frenkel, R. Jack, 
and B. Mladek for discussions at various stages of this 
project. We acknowledge National Science Foundation 
Grant No. NSF DMR-1055586. 



[1] F. Cinti, P. Jain, M. Boninsegni, A. Micheli, P. ZoUer, 

and G. Pupillo, Phys. Rev. Lett. 105, 135301 (2010). 
[2] S. Saccani, S. Moroni, and M. Boninsegni, Phys. Rev. B 

83, 092506 (2011). 
[3] C. N. Likos, Soft Matter 2, 478 (2006). 
[4] D. A. Lenz, B. M. Mladek, C. N. Likos, G. Kahl, and 

R. Blaak, J. Phys. Chem. B 115, 7218 (2011). 
[5] K. Klymko and A. Cacciuto, Phys. Rev. Lett. 107, 

278302 (2011). 
[6] C. N. Likos, A. Lang, M. Watzlawek, and H. Lowen, 

Phys. Rev. E 63, 031206 (2001). 
[7] B. M. Mladek, P. Charbonneau, C. N. Likos, D. Frenkel, 

and G. Kahl, J. Phys.: Condens. Matter 20, 494245 

(2008). 
[8] F. H. Stillinger, J. Chem. Phys. 65, 3968 (1976). 
[9] A. Lang, C. N. Likos, M. Watzlawek, and H. Lowen, J. 

Phys.: Condens. Matter 12, 5087 (2000). 
[10] S. Prestipino, F. Saija, and P. V. Giaquinta, Phys. Rev. 

E 71, 050102 (2005). 
[11] A. Ikeda and K. Miyazaki, Phys. Rev. Lett. 106, 015701 

(2011). 
[12] A. Ikeda and K. Miyazaki, J. Chem. Phys. 135, 024901 

(2011). 
[13] C. Marquest and T. A. Witten, J. Phys.-Paris 50, 1267 



(1989). 
[14] C. N. Likos, M. Watzlawek, and H. Lowen, Phys. Rev. E 

58, 3135 (1998). 
[15] M. Schmidt, J. Phys. Condens Matter 11, 10163 (1999). 
[16] A. Santos, AIP Conf. Proc. 762, 276 (2005). 
[17] G. Falkinger, Master's thesis, Institut fiir Theoretische 

Physik, TU Wien (2007). 
[18] S. H. Suh, C. H. Kim, S. C. Kim, and A. Santos, Phys. 

Rev. E 82, 051202 (2010). 
[19] B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, and 

C. N. Likos, Phys. Rev. Lett. 96 (2006). 
[20] H. Fragner, Phys. Rev. E 75, 061402 (2007). 
[21] B. M. Mladek, P. Charbonneau, and D. Frenkel, Phys. 

Rev. Lett. 99, 235702 (2007). 
[22] A. J. Moreno and C. N. Likos, Phys. Rev. Lett. 99 (2007). 
[23] T. Neuhaus and C. N. Likos, J. Phys.: Condens. Matter 

23 (2011). 
[24] K. Zhang, P. Charbonneau, and B. M. Mladek, Phys. 

Rev. Lett. 105, 245701 (2010). 
[25] A. Nikoubashman, G. Kahl, and C. N. Likos, Phys. Rev. 

Lett. 107, 068302 (2011). 
[26] W. C. Swope and H. C. Andersen, Phys. Rev. A 46, 4539 

(1992). 
[27] The vocable "penetrable sphere model" is also used for 



12 



the Widom-Rowlinson (WR) model, in which overlap- 
ping particles attract with many-body interactions de- 
termined by the total occupied volume [53. 1 54| . The WR 
model, however, thermodynamically maps to a binary 
mixture with non-additive hard-core repulsions [5a|, not 
to the model studied here. [42 

[28] The differential can be obtained for any Uc, but only has [43 
a clear-cut thermodynamic interpretation when its value [44 
is calculated for the equilibrium cluster occupancy. 

[29] R. Kubo, Rep. Prog. Phys. 29, 255 (1966). [45 

[30] M. Allen and D. J. Tildesley, Computer Simulation of [46 
Liquids (Oxford University Press, New York, 1987). 

[31] G. Orkoulas and D. P. Noon, J. Chem. Phys. 131, 161106 [47 
(2009). 

[32] E. A. Guggenheim, J Chem. Phys. 7, 103 (1939). [48 

[33] T. L. Hill, Statistical Mechanics (Dover Publications, 

New York, 1987). [49 

[34] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 

(1992). [50 

[35] A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, 

and P. N. Vorontsovvelyaminov, J. Chem. Phys. 96, 1776 [51 
(1992). 

[36] F. G. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 [52 
(2001). 

[37] F. G. Wang and D. P. Landau, Phys. Rev. E 64 (2001). [53 

[38] D. Frenkel and B. Smit, Understanding Molecular Simu- 
lation (Academic Press, London, 2002). [54 

[39] M. E. J. Newman and G. T. Barkema, Monte Carlo 

Methods in Statistical Physics (Oxford University Press, [55 
New York, 1999). 

[40] A. J. Schultz and D. A. Kofke, Phys. Rev. E 84, 046712 

(2011). [56 

[41] The correct usage of this al gori thm relies on conserving 
the crystal's center of mass [40], for which one needs to 
be particularly careful in cluster crystals [4gl. A possible 



solution would be to track the shift in the center of mass 

using a Fourier transform of the configuration, as was 

done in Ref. ISg . Yet because large, low T cluster crystals 

do not noticeably drift over the simulation time, we here 

neglected this fairly small effect. 

R. D. Kaminsky, J. Chem. Phys. 101, 4986 (1994). 

F. A. Escobedo, J. Chem. Phys. 127, 174104 (2007). 

K. K. Mon and R. B. Griffiths, Phys. Rev. A 31, 956 

(1985). 

J. G. Kirkwood, J. Chem. Phys. 3, 300 (1935). 

B. Mladek, Ph.D. thesis, Institut fiir Theoretische 

Physik, TU Wien (2007). 

W. G. T. Kranendonk and D. Frenkel, Mol. Phys. 72, 

679 (1991). 

P. C. Hemmer, E. Velasco, L. Mederos, G. Navascues, 

and G. Stell, J Chem. Phys. 114, 2268 (2001). 

T. Zykova-Timan, J. Horbach, and K. Binder, J. Chem. 

Phys. 133, 014705 (2010). 

W. G. T. Kranendonk and D. Frenkel, Mol. Phys. 72, 

699 (1991). 

F. Smallenburg, L. Filion, M. Marechal, and M. Dijkstra, 

arXiv: 1111.3466v3 (2012). 

O Frenkel (2011), URL 

|http://www .condmatj ournalclub.org/?p=1619 | 
B. Widom and J. S. Rowhnson, J. Chem. Phys. 52, 1670 
(1970). 

J. S. Rowlinson and B. Widom, Molecular Theory of Cap- 
illarity (Dover Publications, New York, 2003). 
L. Samaj, in Modern Encyclopedia of Mathematical 
Physics: Selecta, edited by I. Y. Aref'eva and D. Stern- 
heimer (Springer, 2007), arXiv: 0709.0617vl. 
K. Zhang and P. Charbonneau, Phys. Rev. B 83, 214303 
(2011). 



